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Abstract 

A  recently-introduced  multiscale  framework  is  used  to  develop  efficient  analysis  and  design  tech-' 
niques  for  networks  with  self-similar  traffic.  These  allow  the  interarrival  density  function  for  fractal 
point  processes  under  Bernoulli  random  erasure  to  be  determined,  as  well  as  the  counting  process 
distribution  for  superpositions  of  these  processes.  The  results  suggest  that  fractal  characteristics 
are  preserved  under  traffic  branching  and  merging,  which  may,  in  turn,  provide  insight  into  the 
prevalence  of  self-similarity  in  aggregate  traffic  broadly  observed  on  real  networks. 

Multiscale  techniques  are  also  developed  for  analyzing  fractal  queueing  scenarios.  These  are  used 
to  obtain,  as  examples,  the  steady-state  customer  distribution  for  a  memoryless  queue  servicing  self¬ 
similar  arrivals,  and  for  Poisson  customers  serviced  with  self-similar  holding  times.  The  persistent 
memory  inherent  in  the  underlying  point  processes  leads  to  substantially  different  behavior  than  is 
observed  in  traditional  queueing  scenarios,  and  important  implications  on  resource  consumption  and 
quality  of  service  are  discussed. 

We  show  how  multiscale  methods  can  be  used  in  conjunction  with  dynamic  programming  tech¬ 
niques  to  develop  efficient  and  practical  control  policies  for  these  fractal  queueing  scenarios.  In 
particular,  optimal  server  control  is  developed  for  a  memoryless  queueing  system  with  self-similar 
traffic  input,  and  optimal  flow  control  is  formulated  for  self-similar  service  of  memoryless  traffic.  By 
exploiting  past  history,  these  controllers  achieve  substantially  better  performance — both  in  terms  of 
quality  of  service  and  resource  utilization — than  traditionally  used  queueing  control  strategies. 

Index  Terms — fractals,  point  processes,  queueing  theory,  optimal  control,  teletraffic,  networks,  self¬ 
similar  traffic,  dynamic  programming 

1  Introduction 


Point  processes  and  queueing  systems  with  fractal  properties  are  increasingly  being  viewed 
as  important  models  in  a  host  of  communication  network  applications.  In  particular,  self¬ 
similarity  in  point  processes  is  well-matched  to  the  burstiness  observed  in  many  aspects  of 
such  networks.  A  brief  listing  includes  error  occurrence  on  a  number  of  telecommunication 
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links,  arrival  patterns  of  many  forms  of  data  such  as  compressed  video,  as  well  as  the  ag¬ 
gregate  traffic  on  a  wide  range  of  networks  like  local  ethernets  and  the  Internet  as  a  whole 
[1,  2,  3,  4].  At  the  same  time,  queueing  models  with  self-similar  properties  are  equally 
promising,  capturing  many  aspects  of  packet  networks  for  which  traditional  memoryless 
models  are  overly  simplistic.  For  example,  power-law  holding  times  of  fractal  queues  are 
well-matched  with  the  propagation  delay  associated  with  heavy-tailed  packet-size  distribu¬ 
tions. 

Motivated  by  the  ubiquity  of  scale-free  event  distributions  in  these  and  a  variety  of  other 
apphcations,  a  number  of  fractal  point  process  models  have  been  developed  and  explored 
in  the  hterature  [3,  5,  6,  7].  In  this  paper,  we  develop  a  multiscale  framework  for  the 
analysis  and  design  of  networks  involving  self-similar  point  processes  and  queues.  Our 
resulting  analysis  results  lead  to  techniques  for  predicting  the  impact  that  the  presence  of 
self-similar  distributions  has  on  the  performance  of  telecommunications  networks.  In  turn, 
our  subsequent  investigation  of  network  design  focuses  on  methods  of  optimal  control  of 
network  activities.  Collectively,  these  results  have  significant  implications  in  the  optimal 
structuring  and  management  of  both  existing  and  future  networks. 

The  rest  of  this  paper  is  organized  as  follows.  In  Section  2,  we  summarize  key  aspects 
of  the  fractal  point  process  model  and  an  associated  multiscale  framework  developed  in  [8], 
which  forms  the  basis  for  our  current  development.  In  Section  3,  we  introduce  a  Marko¬ 
vian  interpretation  for  this  multiscale  framework,  which  is  well  suited  for  determining  the 
counting  process  statistics  for  this  process.  In  Section  4,  we  apply  this  framework  to  the 
analysis  of  key  network  activities  involving  fractal  traffic,  which  include  interaction  among 
multiple  traffic  streams  and  queueing  of  fractal  processes.  In  Section  5,  we  develop  exten¬ 
sions  of  our  multiscale  framework  to  allow  us  to  design  control  policies  for  a  number  of 
fractal  queues,  including  optimal  server  and  flow  control.  Finally,  Section  6  contains  some 
concluding  remarks. 


2  Fractal  Renewal  Processes 


A  point  process — i.e.,  a  random  distribution  of  event  arrivals  in  time — is  naturally  chax- 
acterized  in  terms  of  its  interarrival  intervals  A’[i].  More  precisely,  we  let  X[l]  denote  the 
arrival  epoch  of  the  first  arrival  after  t  =  0,  and  A'[i]  denote  the  time  interval  between  the 
{i  —  l)st  arrival  and  the  ith  arrival,  for  i>2.  Other  aspects  of  a  point  process  axe  revealed 
by  its  characterization  in  terms  of  the  associated  counting  process  N{t),  whose  value  at 
time  t  is  defined  as  the  number  of  arrivals  in  the  interval  (0,  t] .  Since  the  counting  process 
describes  the  history  of  the  point  process,  its  value  at  any  instant  t  has  dependence  on  the 
choice  of  the  reference  point  t  =  0.  Two  scenarios  are  of  special  interest.  In  the  arrival- 
observed  case,  the  reference  point  coincides  with  an  event.  In  the  context  of  networking 
studies,  for  example,  this  can  represent  data  traffic  as  viewed  by  a  user.  On  the  other  hand, 
the  random  incidence  case  corresponds  to  the  point  of  reference  being  chosen  randomly, 
independent  of  the  point  process.  As  such,  random  incidence  is  useful  for  modeling  traffic 
firom  the  perspective  of  network  management,  for  example. 

The  fractal  point  process  model  of  interest  is  defined  in  terms  of  a  particular  self- 

'p 

similarity  property.  Specifically,  the  associated  coimting  process  satisfies  N{t)  =  N{at)  for 

•p 

all  o  >  0,  where  the  notation  =  denotes  statistical  equality,  in  particular  in  the  sense  of  all 
finite-dimensional  distributions. 

Much  of  the  physical  network  behavior  of  interest  is  effectively  stationary,  exhibiting 
no  preference  for  a  time  origin.  As  such,  it  is  tempting  to  restrict  attention  to  self-similar 
point  processes  that  are  also  renewal  processes.  However,  a  true  renewal  process  cannot  be 
self-similar  [8].  Nevertheless,  it  is  possible  to  develop  a  fractal  point  process  model  based 
on  a  generalized  notion  of  renewal  process.  To  develop  this  notion,  we  first  introduce  the 
following  convenient  terminology;  we  say  that  a  point  process  with  interarrivals  y[i]  is 
derived  from  a  point  process  with  interarrivals  X[i]  via  conditioning  on  the  event  S  if  F[i] 
is  the  subsequence  of  X[z]  formed  by  discarding  those  components  such  that  A’[i]  ^  £.  As 
developed  in  [8],  a  fractal  renewal  process  is  then  defined  as  a  self-similar  point  process  that 
satisfies  the  following: 

1.  When  conditioned  on  the  event  £  =  {x  <  JA  <  x}  for  some  0  <  x  <  x  <  oo,  the 
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resulting  point  process  is  a  renewal  process;  and 


2.  When  conditioned  on  each  of  any  number  of  arbitrary,  mutually  exclusive  events 
^1,^2-,  ■■■  iJ-L  such  that  J-i  =  {^  <  X  <  xi)  with  0  <  ^  <  x;  <  oo,  the  resulting 
point  processes  axe  mutually  independent. 

Since  observations  of  physical  signals  are  typically  limited  both  in  resolution  and  duration, 
Property  1  implies  that  a  fractal  renewal  process  is  effectively  a  renewal  process.  Further¬ 
more,  it  can  be  shown  that  as  a  consequence  of  self-similarity,  the  interarrivals  y[i]  resulting 
from  such  observations  must  be  distributed  as  [8] 


/y(y)  x<y<x 


(1) 


where  x  and  x  are  determined  by  the  resolution  and  duration  windowing,  respectively.  The 
shape  parameter  7  in  (1)  is  directly  related  to  the  fractal  dimension  D  of  the  process  via 
7  =  D  -t- 1.  Fig.  1  illustrates  a  typical  sample  function  of  a  fractal  renewal  process  with 
shape  parameter  7  =  1.5,  viewed  under  successive  magnification,  from  which  the  hallmark 
scale-independent  clustering  behavior  is  apparent. 

Analysis  of  the  fractal  renewal  process  is  facilitated  by  a  highly  efficient  Poisson-based 
multiscale  representation  developed  in  [8].^  The  essence  of  this  framework  is  to  model 
a  fractal  renewal  process  as  a  mixture  of  a  multiscale  family  of  Poisson  processes.  In 
terms  of  interarrival  statistics,  this  is  equivalent  to  decomposing  the  power-law  interarrival 
distribution  of  (1)  as  a  weighted  sum  of  dilated  and  compressed  exponential  functions.  In 
particular,  with  a  finite  number  of  constituents,  the  interarrival  decomposition  takes  the 
form 

i-i 

fxix)  =  '^PjfXjix),  (2) 

j=0 

^In  fact,  such  representations  naturally  lead  to  computationally-efEcient  cJgorithms  for  robust  estimation 
of  the  scaling  parameters  of  such  processes,  as  also  developed  in  [8]. 


where  the  Xj  are  exponential  random  variables  with  rate  parameters 


while  the  pj  form  a  geometric  probability  distribution 

Pj  =  >  (4) 

with  cr^  as  a  normalization  constant.  Spacing  between  the  constituents  is  governed  by  the 
scale  increment  rj>  1,  while  the  number  of  scales  is  equal  to  L,  with  scale  0  being  the  finest 
scale. 

3  The  Multiscale  Pure-Birth  Process 

The  preceding  multiscale  model  has  a  natural  interpretation  in  the  form  of  a  multiscale  pure- 
birth  process.  This  process,  depicted  in  Fig.  2,  forms  the  basis  of  our  development.  General¬ 
izing  the  well-known  pure-birth  process  (see,  e.g.,  [9]),  the  state  space  of  this  Markov  process 
is  naturally  partitioned  into  “superstates,”  each  of  which  represents  a  certain  number  of 
births.  A  superstate  is,  in  turn,  composed  of  a  set  of  states  which  correspond  to  the  scales 
in  our  finite-scale  framework.  Hence,  we  index  the  states  with  a  pair  of  integers  {i,j),  where 
the  superstate  index  i  is  nonnegative,  while  the  scale  index  j  ranges  from  0  to  i  —  1  for  an  L- 
scale  representation.  We  define  the  probability  vector  Pj (t)  =  Pi^i(t), . . . ,  Pi^L^i{t)), 

with  Pij{t)  denoting  the  probability  that  the  process  is  in  state  {i,j)  at  time  t.  With  this 
notation,  the  counting  process  probability  distribution  can  be  expressed  as 

Pr{iV(t)  =pi(t)l'^,  (5) 

where  1  is  a  row  vector  with  all  entries  equal  to  1.  For  convenience,  we  use  the  notation 
■ni{t)  =  Pr{Ar(t)  =  i}. 

Every  transition  in  the  multiscale  pure-birth  process  results  in  an  increment  in  the 
superstate  index,  and  thus  a  birth.  In  our  model,  the  birth  rate  is  assumed  independent 
of  the  number  of  births  already  occurred.  Thus,  the  mean  departure  rate  from  each  state 
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(i,j)  is  only  a  function  of  the  scale  j,  taking  the  specific  form  Xj  =  X/rf  in  terms  of  the 
scale  increment  77  and  the  reference  rate  A.  Upon  a  birth,  every  state  f  of  the  succeeding 
superstate  can  be  immediately  reached,  with  probability  pj/  =  cr^q^' ,  where  q  =  governs 
the  relative  weighting  on  the  scales,  while  <7^  is  the  normalization  factor.  Thus,  the  duration 
between  consecutive  births  is  governed  by  a  probability  law  of  the  form  (2). 

The  counting  process  associated  with  the  firactal  renewal  process  can  be  readily  obtained 
from  the  multiscale  pure-birth  representation.  Inspecting  Fig.  2,  it  is  straightforward  to  set 
up  the  forward  Kolmogorov  equations  [9]  for  the  multiscale  pure-birth  process.  In  vector 
form,  we  have  that 

=  -po(t)B  (6a) 

=  -pi(t)B  ■+-pi_i(<)b'^q,  z>l  (6b) 

where  q  =  {po,pi,  ■  ■  ■  ,Pl-i)  contains  the  choice  probabilities,  b  =  is 

a  vector  of  the  multiscale  dilation  factors,  and  B  =  diag(b)  is  a  diagonal  matrix  with  the 
dilation  factors  along  its  main  diagonal.  Prom  (6)  we  obtain  the  transform-domain  equation 

\  ^)  =  P(2;  t)  (-B  -h  zb'^q) ,  (7) 

where  p{z;  t)  denotes  the  z-transform  of  the  sequence  {pi(t);  i  =  0, 1, . . .  },  defined  as 

00 

p(^u)  = 

i=0 

Eq.  (7)  can  now  be  readily  solved  to  yield 

p(2r;  t)  =  p{z;  0)  exp  ([— B  -I-  2:b’^q]  At)  ,  t>  0,  (8) 

which  determines  the  z-transform  up  to  an  initial  condition. 

In  the  arrival-observed  case,  where  the  point  of  reference  t  =  0  is  a  renewal,  the  first 
interaxrival  is  statistically  identical  to  every  other  interarrival.  Hence,  the  scales  within 
the  zeroth  superstate  are  chosen  with  probabilities  {pj;j  =  0, 1, . . . ,  L  —  1},  and  the  initial 


Id  ,  , 
A*""'*’ 

A*"*'*) 


condition  in  (8)  is  p(2;0)  =  q.  Inverting  the  resulting  transform  p(2:;f)l'^,  we  can  obtain 
a  time-domain  characterization  of  the  arrival-observed  counting  process  distribution,  which 
we  denote  by  | (t);  z  =  0, 1, . . .  In  particular, 


L-l 

=  Po(*)l'^  =  P(0;i)l'^  =  qexp(-BAt)l'^  =  Y^piexp(-Xit) ,  (9) 

i=0 


which  decays  as  for  large  t  when  L  oo.  Higher-order  terms  can  be  obtained 

numerically  from  a  Taylor  series  expansion  of  the  transform  p(2r;t)l^,  as  developed  in 
detail  in  Appendix  A.  The  main  results  are  that  for  i  >  1, 


i  “i+i  '  • 


i+1 


{i  -f  1)! 


-h' 


with  the  first-order  coefficients  given  by 


k 

=^J2MiMk-i  (10) 

z=i 

and  the  higher-order  coefficients  obtained  via  the  recurrence  relation 

fc— i-i-i 

0?  =  E  (11) 

1=1 

where  Mi  is  the  1th  moment  of  a  random  variable  R  distributed  according  to 


PT{R  =  ri  ^}=pj,  j  =  0,1,..., L-l. 


(12) 


For  random  incidence,  we  assume  observation  begins  at  a  random  time,  with  the  point 
process  already  in  equilibrium.  Under  this  assumption,  the  scale  of  the  first  interarrival  is 
selected  with  the  steady-state  marginal  distribution  over  scales,  which  is  d^qB~^,  where 
is  for  normalization;  an  argument  is  given  in  Appendix  B.  Using  this  in  (8),  we  can  obtain  a 
time-domain  characterization  of  the  random  incidence  counting  process  distribution  which 
we  denote  by  |7r|'^^(t);i  =  0, 1, . . .  |.  For  example,  it  is  straightforward  to  get  the  closed- 
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form  expression 


i=0 


(13) 


which  when  L  oo  approaches  1  if  7  <  2  and  decays  like  l/f^  ^  for  large  t  if  7  >  2. 
Higher-order  terms  of  the  probability  distribution  take  the  form 


_  (0  , 


where  the  coefficients  are  obtained  from  those  in  the  arrival-observed  distribution  via 


ri*^  =  (4-i+4' 


(14) 


Details  of  the  derivation  axe  given  in  Appendix  B. 

Figs.  3  and  4  illustrate  the  time  evolution  of  the  lower  order  terms  in  the  arrival-observed 
and  random  incidence  counting  process  distributions  for  the  case  7  =  1.8.  It  is  worth 
remarking  that  in  practice,  very  few  scales  are  typically  needed  for  a  good  approximation 
to  the  probabilities  over  a  finite  time  interval,  although  more  scales  are  generally  required 
in  small  7  situations  due  to  the  more  persistent  tail. 

Several  features  of  the  plots  are  noteworthy.  First,  in  the  arrival-observed  case  of 
Fig.  3,  that  the  probabilities  Tri{t)  all  fall  oflF  along  the  same  asymptote  implies 

that  £J[A^(t)]  ~  These  statistics  axe  consistent  with  the  strong  clustering  behav¬ 

ior  characteristic  of  these  processes,  which  is  increasingly  pronounced  as  7  decreases.  By 
contrast,  for  the  memoryless  (Poisson)  process,  the  counting  process  probabilities  fall  ofiT 
exponentially  quickly  and  E  [A^(t)]  ~  0{t).  Moreover,  the  power-law  probability  decay  is 
also  consistent  with  the  results  from  direct  computation  of  the  counting  process  statistics 
[10],  which  involves  asymptotic  analysis  of  convolution  of  multiple  power- law  functions.  In 
addition,  the  counting  process  distribution  depicted  in  Fig.  4  for  the  random  incidence  case 
provides  dramatic  evidence  of  the  impact  of  the  unusually  long  quiescent  periods  between 
clusters  in  fractal  point  processes. 

We  turn  now  from  the  analysis  of  firactal  point  processes  in  isolation  to  analysis  of  their 
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behavior  in  scenarios  representative  of  those  encountered  by  traffic  in  large,  interconnected 
networks. 

4  Network  Analysis  for  Fractal  Traffic 

Key  activities  in  a  communication  network  can  typically  be  modeled  in  terms  of  transfor¬ 
mations  on  point  processes.  Interaction  among  traffic  streams,  for  example,  can  often  be 
modeled  as  branching  and  merging  of  point  processes.  Propagation  delay  and  processing 
latency  at  sites  and  hnks  can  typically  be  captured  by  queueing  models.  In  this  section, 
we  apply  multiscale  techniques  to  analyze  a  number  of  these  activities.  Branching  and 
erasure  of  fractal  renewal  processes  will  be  addressed  in  Section  4.1,  while  fractal  renewal 
process  superposition  will  be  the  theme  of  Section  4.2.  Queueing  systems  with  self-similar 
properties  will  be  explored  in  Sections  4.3  and  4.4. 

4.1  Random  Erasure  of  the  Fractal  Renewal  Process 

Data  loss  and  branching  of  data  streams  are  two  of  the  most  frequently  encountered  activ¬ 
ities  in  a  network.  An  often  reahstic  and  widely-adopted  model  for  both  transformations 
is  Bernoulli  point  process  erasure,  whereby  each  point  is  independently  erased  with  a  com¬ 
mon  probability  p.  In  what  follows,  we  study  the  behavior  of  the  fractal  renewal  process 
under  this  mode  of  erasure.  A  key  result  of  our  analysis  is  the  preservation  of  self-similar 
characteristics  even  under  very  high  erasure  probabilities. 

To  determine  interarrival  density  under  Bernoulli  erasure,  we  begin  by  observing  that 
with  probability  1  —  p,  an  arrival  of  the  original  process  contributes  a  count  of  unity  to 
the  counting  process  resulting  from  erasure.  Otherwise,  it  contributes  a  count  of  zero.  In 
the  transform  domain,  this  corresponds  to  the  usual  replacing  of  z  with  p  -h  (1  -  p)z  in 
the  transform  of  the  original  counting  process  distribution.  In  particular,  using  (8),  the 
counting  process  distribution  of  an  erased  fractal  renewal  process  has  the  z-transform 

p(p-t-  (1  -p)z;t)i^  =  p(2;;0)exp  ([-B  -I-  (p-l-  (1  -p)z)b'^q]  At)  l'^.  (15) 

With  the  initial  condition  p(2:;0)  =  q,  and  with  =  0  in  (15),  the  arrival-observed  proba- 
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bility  of  zero  arrivals  in  an  interval  (0,  w]  is 


7ro®'^(w)  =  qexp  ([-B  +  pb'^q]  Aiu)  l'^.  (16) 

But  this  event  is  equivalent  to  the  event  {W  >  i/;},  where  W  is  the  interarrival  beginning 
at  0.  DiflFerentiating  (16),  we  have  the  interarrival  density 

fwiw)  =  Pr{W'  >w}  =  qexp  ([— B  +pb^q]  Xw)  (B  —  pb"^q) 

=  A(1  —  p)qexp  ([— B +pb^q]  Atu)  b"^,  to  >  0.  (17) 

Using  (17),  we  have  plotted  in  Fig.  5  the  interarrival  density  of  a  fractal  renewal  pro¬ 
cess  with  shape  parameter  7  =  1.8,  subject  to  various  erasure  probabilities.  These  plots 
suggest  that  the  erased  interarrival  density  largely  retains  the  power-law  characteristics  of 
the  original  density  (i.e.,  the  p  =  0  case).  In  fact,  for  p  <  1,  empirical  studies  suggest  that 
as  A  — >•  00,  L  00, 

Mw)  ~  i 


for  every  lu  >  0. 

4.2  Superposition  of  Fractal  Renewal  Processes 

Merging  of  data  streams  is  another  point  process  transformation  typical  in  networks.  To 
investigate  the  behavior  of  fractal  point  processes  under  merging,  we  consider  the  super¬ 
position  of  two  independent  fractal  renewal  processes.  A  key  implication  of  our  results  is 
the  invariance  of  fractal  point  process  features  under  superposition.  More  importantly,  our 
results  also  suggest  that  the  family  of  fractal  point  processes  constitutes  a  domain  of  attrac¬ 
tion  under  aggregation,  much  like  the  Poisson  family.  This  behavior  is  consistent  with  the 
spectral  analysis  results  obtained  in  [6].  Together  with  the  random  erasure  results,  these 
superposition  results  may  prove  useful  in  explaining  empirical  observations  of  the  ubiquity 
of  self-similarity  in  aggregate  traffic  on  a  broad  range  of  networks. 

When  two  point  processes  are  superimposed,  their  counting  processes  add.  Prom  an 
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arrival’s  viewpoint,  the  structure  of  the  two  counting  process  distributions  are  inherently 

different.  On  one  hand,  the  process  to  which  it  belongs  is  governed  by  the  arrival-observed 

distribution;  we  denote  this  distribution  by  i  =  0, 1, . . .  On  the  other  hand, 

since  the  constituents  are  independent,  the  arrival  observes  the  random  incidence  counting 

processing  distribution,  denoted  as  i  =  0, 1, . . .  j-,  for  the  other  point  process. 

fl  a) 

Thus,  the  overall  counting  process  distribution  will  be  a  discrete  convolution  of  ’  ’{t) 
and  If,  in  addition,  the  two  constituents  have  the  same  fractal  dimension,  the 

resulting  counting  process  distribution  will  be  expressed  more  simply  as 

W  =  ^  i  =  0, 1, . . . .  (18) 

j=0 

By  a  similar  argument,  the  random  incidence  counting  process  distribution  of  the  superpo¬ 
sition  is  the  convolution  of  the  two  constituent  random  incidence  counting  process  distri¬ 
butions, 

*  =  0, 1, ... .  (19) 

j=o 

Figs.  6  and  7  show  the  arrival-observed  and  random  incidence  counting  process  distribu¬ 
tions  corresponding  to  the  superposition  of  two  independent  fractal  renewal  processes,  each 
with  shape  parameter  7  =  1.8.  The  computations  were  performed  according  to  (18)  and 
(19),  respectively,  using  the  counting  process  distribution  results  of  Section  3.  Comparing 
this  set  of  plots  with  those  for  a  single  process  (Figs.  3  and  4),  we  observe  that  key  features 
such  as  the  asymptotic  power-law  decay  in  the  arrival-observed  distribution,  and  dominance 
of  zeroth-order  term  in  the  random  incidence  distribution,  are  largely  preserved  under  su¬ 
perposition,  suggesting  invariance  of  fractal  renewal  processes  under  this  transformation. 
These  counting  process  results  provide  additional  verification  of  the  invariance  of  firactal 
point  proceses  under  superposition,  complementing  the  spectral  evidence  given  in  [6] . 
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4.3  Queueing  of  the  Fractal  Renewal  Process 

While  random  erasme  and  superposition  capture  key  interaction  among  multiple  data 
streams,  queueing  models  are  appropriate  for  activities  at  individual  sites  and  channels. 
More  generally,  queueing  analysis  is  important  in  many  applications  involving  resource 
sharing.  The  multiscale  pure-birth  model  can  be  readily  extended  to  support  queueing  sce¬ 
narios.  For  example,  addition  of  backward  transitions,  or  deaths,  results  in  the  multiscale 
birth-and-death  process  of  Fig.  8,  which  models  a  fractal  renewal  process  being  serviced 
by  a  single-server,  memoryless  queue.  The  added  transitions  model  service  completion  and 
customer  departure  and  therefore  all  occur  at  the  mean  service  rate  /i.  Because  the  ser¬ 
vice  being  modeled  is  memoryless,  the  multiscale  birth-and-death  model  is  applicable  to 
work-conservative  queues  with  arbitrary  service  discipline. 

As  an  illustration  of  the  types  of  results  that  can  be  obtained  with  this  Markov  model, 
we  develop  systematic  and  efficient  methods  to  compute  the  steady-state  probability  dis¬ 
tribution  of  the  number  of  customers  in  the  system  (tti;  i  =  0, 1, . . .  },  which  is  particularly 
insightful  for  buffer  allocation  and  predicting  quality  of  service,  for  example.  The  dynamics 
of  the  process  in  Fig.  8  are  governed  by  the  system  of  forward  Kolmogorov  equations: 

~po(t)  =  -po(t)B  +  -pi(t)I  (20a) 

Xat  p 

“Pi(<)  =  -f  -I-  pi_i(i)b'^q  +  ipi+i(t)I,  i  >  1,  (20b) 

where  p  =  X/p.  To  obtain  the  equilibrium  state  distribution  (pj;  z  =  0, 1, . . .  },  we  employ 
the  matrix-geometric  methods  of  [11].  By  this,  we  first  obtain  a  positive  semi-definite 
solution  R  to  the  quadratic  equation 

0  =  pb^q  -  R(pB  -I- 1)  -I-  R2.  (21) 

This  can  be  accomplished  via  successive  approximation  (see  [11]).  The  steady-state  proba¬ 
bilities  are  then  given  by 

Pi  =  xR*,  i  =  0, 1, . . . ,  (22) 
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where  x  is  the  (unique)  left  null  vector  of  pB  —  R,  normalized  by  the  constraint 

x(I  -  =  1.  (23) 

Finally,  the  steady-state  customer  distribution  can  be  obtained  via  =  pjl^. 

These  methods  axe  a  much  more  computationally  practical  means  for  predicting  queue 
performance  than  are  simulation-based  approaches.  In  Fig.  9,  the  analytical  method  for 
determining  the  steady-state  distribution  for  the  nmnber  of  customers  in  the  system  with 
service  rate  l/p  =  p/A  =  0.15  and  shape  parameters  7  =  1.5, 1.8  is  compared  with  that 
obtained  by  simulations  using  the  next-event  time  advance  methods  of  [12].  Despite  large 
sample  sizes,  discrepancy  between  the  theoretical  and  simulation  results  remains  largely  as 
a  consequence  of  substantial  outliers  in  the  simulations.  Such  outliers  are  characteristic  of 
the  heavy-tailed  distributions  involved. 

Both  the  theoretical  and  simulation  results  suggest  that  conditioned  on  one  or  more 
customers,  the  distribution  approaches  a  geometric  function.  This  is  indeed  the  case,  and 
can  be  attributed  to  the  fact  that  the  solution  to  (21)  is  of  rank  1.  In  particular,  it  can  be 
shown  that  the  solution  must  be  of  the  form 


R  =  b^w. 


(24) 


which  implies  that  the  steady-state  customer  distribution  satisfies 

=  wb^  =  A:,  t  =  1, 2, . . . .  (25) 

TTi 

Thus,  k  is  the  decay  rate  of  the  geometric  portion  of  the  distribution.  From  (21),  we  get 
that  w  in  (24)  must  satisfy 

pq  =  w  (pB  -I- 1)  —  wfe  =  w  [pB  -f  (1  —  k)  I] . 
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Right  multiplication  of  both  sides  by  [pB  +  (1  —  A:)I]  ^  leads  to 


+ 


Qph 


r,L-l 


lp+{l-k)  p/r)  +  {l-k) 


p/t] 


L-l 


p/ri^~^  +  (1  -  fe) 


=  k. 


(26) 


which  can  be  solved  iteratively  using  bisection  search  techniques,  for  instance. 

A  useful  bound  for  k  can  be  obtained  via  Jensen’s  inequality.  Specifically,  the  left  hand 
side  of  (26)  is  less  than  or  equal  to 


pE[R] 

pE[R]  +  {i-ky 


(27) 


where  R  is  the  random  variable  defined  in  (12).  Simplifying,  we  get  that 


A:(l-fc)  <  pE[R](l-k), 


which,  since  A:  <  1  for  an  ergodic  queue,  simplifies  to 

k<pE[R]  =  ^E[R].  (28) 

Fig.  10  shows  values  of  k  obtained  via  bisection-search  solution  of  (26),  together  with  the 
bounds  as  prescribed  by  (28).  As  expected,  higher  service  efficiency  (p/A)  leads  to  less  con¬ 
gested  system,  as  reflected  by  the  lower  k  (or,  sharper  decay  in  the  customer  distribution). 
These  plots  also  suggest  that  while  the  boxmd  in  (28)  is  somewhat  loose  for  slow-service 
scenarios,  it  is  a  potentially  useful  closed-form  approximation  for  k  in  the  fast-service  regime. 

Another  important  feature  of  the  distributions  in  Fig.  9  is  the  dominance  of  the  prob¬ 
ability  of  zero  customers  (i.e.,  idle  system).  This  can  again  be  attributed  to  the  unusually 
long  gaps  between  clusters  of  customer  arrivals.  This  behavior  suggests  that  service  rate 
may  be  lowered  for  large  time  intervals,  with  no  noticeable  degradation  in  the  quality  of 
service.  In  Section  5.1,  we  show  that  this  can  indeed  be  accomplished  via  dynamic  server 
control. 


4.4  Power-Law  Service  of  the  Poisson  Process 


While  the  system  of  Section  4.3  addresses  self-similarity  in  the  arrivals,  fractal  characteristics 
can  often  be  found  in  other  aspects  of  a  queueing  scenario.  For  example,  the  distribution 
of  packet  sizes  in  many  networks  is  often  well-modeled  as  a  power  law,  with  frequent  short 
packets  and  occasional  extraordinarily  long  packets.  The  power-law  holding  time  arising 
for  this  packet  distribution  can  be  modeled  with  the  Markov  process  of  Fig.  11,  obtained  by 
transposing  the  multiscale  birth-and-death  process  of  Fig.  8.  Thus,  the  input  process  is  now 
Poisson  with  arrival  rate  A,  while  the  service  is  captured  by  an  L-scale  representation  with 
finest-scale  service  rate  fi.  Since  the  service  is  no  longer  memory  less,  this  model  is  restricted 
to  queueing  disciplines  with  no  preemption.  As  such,  various  forms  of  time-slicing  are 
precluded,  for  example.  Nevertheless,  many  important  disciplines  axe  still  captmred,  such 
as  first-in-first-out,  last-in-first-out,  and  service  in  random  order. 

To  obtain  the  steady-state  customer  distribution  of  this  queueing  system,  we  again  use 
the  matrix-geometric  method  of  [11].  In  particular,  we  solve  the  matrix  quadratic  equation 

0  =  pi  -  R(pl  +  B)  -f-  R^b^q  (29) 

for  R  via  successive  approximation,  where  as  before,  p  =  A/p.  Once  R  is  obtained,  the 
steady-state  customer  distribution  can  be  generated  via  (22)  and  (23).  However,  unlike 
(21)  in  the  system  of  Section  4.3,  the  solution  to  (29)  is  not  rank  1,  precluding  further 
simplication. 

Two  sets  of  results  are  plotted  in  Fig.  12,  based  on  the  above  computaton  and  next- 
event  time  advance  simulation.  For  the  top  curves,  the  input  arrival  rate  is  p  =  2  x  10“^, 
while  the  shape  parameter  of  the  service  duration  is  7  =  1.5.  For  the  bottom  curves,  the 
arrival  rate  is  set  at  p  =  1  x  10~®,  while  the  shape  parameter  of  the  service  time  is  7  =  1.8. 
Agreement  between  our  prediction  and  simulation  results  is  close. 

The  figure  also  suggests  that  a  key  feature  of  this  queueing  scenario  is  the  heavy-tailed 
customer  distribution,  which  asymptotically  approaches  a  power  law.  An  imphcation  of 
this  is  the  requirement  of  large  buffers  for  its  implementation.  In  addition,  long  delay  and 
hence  poor  quality  of  service  can  be  expected  for  the  customers  if  a  first-in-first-out  (FIFO) 


discipline  is  adopted.  In  Section  5.2.1,  we  develop  an  optimal  dynamic  queueing  control 
strategy  to  improve  performance  of  this  system. 

5  Network  Control  for  Fractal  Traffic 

The  analysis  in  Sections  4.3  and  4.4  focusses  on  the  behavior  of  queues  with  fixed  system 
parameters.  However,  constant-rate  memoryless  service  of  input  processes  with  self-similar 
arrivals  can  result  in  tremendously  inefiicient  utilization  of  resources,  as  reflected  by  very 
high  idle  probabihty.  Similarly,  poor  quality  of  service  can  result  for  service  of  a  constant- 
rate  Poisson  input  when  the  holding-time  statistics  are  self-similar,  which  manifests  itself 
in  the  form  of  a  heavy-tailed  queue-length  distribution. 

In  many  realistic  scenarios,  however,  various  aspects  of  a  queueing  system  are  control¬ 
lable,  often  in  real  time.  For  example,  controllability  of  service  rate  is  quite  feasible.  In 
communications  engineering,  a  number  of  schemes  have  been  proposed  for  versatile  alloca¬ 
tion  of  bandwidth,  ranging  from  fast  packet  switching  networks  [13]  to  flexible  assignment 
of  multilevel  trunks  and  trunk  groups  [14].  Likewise,  control  of  input  arrival  rate,  or  flow 
control,  is  also  possible  in  many  situations.  Specifically,  input  throttling,  or  arrival  rate 
reduction,  can  be  implemented  via  admission  toll  or  traffic  re-routing. 

Dynamic  programming  [15, 16]  is  a  natural  tool  for  developing  queueing  policies.  Indeed, 
dynamic  programming  has  been  used  to  develop  a  number  of  optimal  controllers  for  the 
M/M/m  family  of  queues,  i.e.,  queueing  systems  with  memoryless  input  and  finitely  many 
(m)  memoryless  servers.  For  example,  it  is  straightforward  to  develop  server  control  for 
the  M/M/1  queue  that  optimizes  a  long-term  cost  of  two  components,  corresponding  to  the 
costs  of  service  and  buffer  occupancy,  respectively  [17].  In  this  scenario,  it  is  well-known  that 
the  optimal  service  rate  at  any  time  depends  only  on  the  number  of  customers  in  system. 
Further,  under  rather  broad  conditions,  the  optimal  service  rate  increases  monotonically 
with  the  number  of  customers  present. 

Similar  techniques  can  be  used  to  design  optimal  flow  control  for  the  M/M/1  queue  [17]. 
In  this  case,  the  admission  rate  is  designed  to  minimize  a  cost  function  composed  of  a  cost 
on  buffer  occupancy  and  a  penalty  for  input  throttling.  As  in  the  server  control  problem. 


optimal  strategy  for  this  problem  depends  only  on  the  number  of  customers  present,  and 
the  results  agree  closely  with  intuition:  imder  rather  broad  conditions,  customer  arrivals 
are  increasingly  deterred  as  the  queue  becomes  more  crowded. 

That  the  optimal  server  and  flow  control  policies  for  M/M/m  systems  do  not  exploit 
history  of  the  system  and  are  completely  determined  by  the  current  state  is  a  direct  con¬ 
sequence  of  the  memoryless  natme  of  such  systems.  In  contrast,  we  show  in  this  section 
that  optimal  control  strategies  for  fractal  queues  depend  heavily  on  past  history  to  main¬ 
tain  efficient  operation  due  to  the  long-term  dependence  in  self-similar  point  processes. 
These  controllers,  which  we  develop  by  applying  dynamic  programming  techniques  within 
our  multiscale  modeling  framework,  significantly  improve  performance  in  terms  of  both 
resource  usage  and  quality  of  service  over  controllers  that  do  not  exploit  past  history. 

5.1  Server  Control  for  the  Fractal  Renewal  Process 

We  first  consider  the  optimal  control  of  the  queueing  system  of  Section  4.3,  where  a  single 
memoryless  server  processes  a  fractal  renewal  process  input.  As  is  common  in  other  state- 
space  control  problems,  we  develop  optimal  server  control  policies  for  this  system  based  on  a 
two-step  approach;  we  first  develop  ideal  controllers  that  rely  on  complete  state  information, 
then  develop  practical  controllers  by  replacing  the  state  information  with  suitably  defined 
state  estimates.  The  first  of  these  subproblems  is  addressed  in  Section  5.1.1;  the  second  is 
treated  in  Section  5.1.2. 

5.1.1  State-Based  Multiscale  Server  Control 

To  model  the  behavior  of  this  system,  we  again  employ  the  multiscale  birth-and-death  model 
of  Fig.  8.  However,  the  service  rate  (j,  is  now  assumed  controllable  and  can  be  varied  over  an 
achievable  range  0  <  fj,  <Ji  based  on  the  state  of  the  queue.  To  reflect  its  state  dependence, 
we  use  the  notation  fuj  to  denote  the  service  rate  when  i  customers  are  present  and  the 
scale  of  the  next  interaxrival  is  j  [i.e.,  when  the  state  is  (i,  j)]. 

A  fundamental  tradeoff  exists  in  the  operation  of  this  queue.  On  one  hand,  as  pointed 
out  in  Section  4.3,  high  service  rates  will  result  in  inefficiently  used  server,  as  reflected  by  the 
substantial  idle  probability.  On  the  other  hand,  if  the  service  rate  is  excessively  reduced. 
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the  quality  of  service  can  severely  decline.  In  the  extreme  case,  the  system  will  become 
non-ergodic  and  the  queue  will  grow  without  bound. 

To  facilitate  quantitative  analysis,  we  assign  costs  to  these  two  criteria.  As  a  measure  of 
resomrce  consumption,  we  define  a  service  cost  c(/z),  which  is  a  continuous,  nondecreasing 
function  of  the  service  rate.  For  convenience,  an  idle  server  is  assumed  to  inflict  zero 
cost,  i.e.,  c(0)  =  0.  To  capture  quality  of  service,  we  define  a  holding  cost  h{i),  which  is 
nondecreasing  in  i,  the  number  of  customers  in  the  system.  This  cost  is  directly  related  to 
the  waiting  time  experienced  by  each  customer  under  the  first-in-first-out  (FIFO)  discipline, 
for  example.  For  convenience,  an  empty  system  is  assumed  cost-firee,  i.e.,  /i(0)  =  0.  The 
overall  objective  function  is  then  the  expected  value  of  the  combined  costs,  accumulated 
over  time: 


J  =  E 


■  poo 

Ja  J 


(30) 


where  i{t)  is  the  number  of  customers  in  the  queueing  system  at  time  t,  and  fi{t)  is  the 
service  rate  at  time  t.  The  discount  rate  (3  is  included  to  allow  weighting  of  future  costs 
with  respect  to  the  present. 

To  exploit  results  firom  discrete-time  dynamic  programming,  we  recast  our  continuous¬ 
time  Markov  decision  problem  into  an  equivalent  discrete-time  one.  We  begin  by  adding 
self-transitions  in  the  multiscale  birth-and-death  process  so  that  the  total  departure  rate 
firom  each  state  is  =  A  H-p,  i.e.,  the  maximum  transition  rate  in  the  original  process.  This 
yields  the  equivalent  Markov  process  shown  in  Fig.  13.  Since  the  departme  rates  axe  now 
uniform  across  all  states,  this  step  is  sometimes  referred  to  as  rate  uniformization  [17]. 

Next,  we  rewrite  the  cost  in  (30)  as 


J  =  E 


ftn+l 

Y;  e-»<it(Mi[n])+c(MH)) 

.n=0 


(31) 


where  {tn',n  =  0, 1, . . .  }  are  the  transition  epochs  of  the  uniformized  process,  and  i{n]  and 
/i[n]  are  respectively  the  number  of  customers  and  the  service  rate  during  the  nth  interval, 
(t„,tn+i].  Carrying  out  the  integration  in  (31)  and  taking  expectation,  we  get  that  the  cost 
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IS 


J=T,E 


n=0 


o  ft  n 


1  °° 


/3 


E[h{i[n])  +c{^J.[n])] 


-E  [h{i[n])  +  c(M[n])] , 


(32) 


where  the  second  equality  follows  from  the  fact  that  tn  is  an  nth-order  Erlang  random 
variable  with  mean  n/Q..  The  final  expression  in  (32)  is  the  objective  function  of  a  discrete¬ 
time  dynamic  programming  problem  for  a  Markov  chain  with  the  same  topology  as  the 
process  in  Fig.  13,  where  for  any  pair  of  states  x  and  y,  the  transition  rate  tx,y  is  replaced 
by  transition  probability  Px,y  =  tx,y/^-  The  holding  and  service  rate  costs  are  respectively 
h{-)l{P  +  f2)  and  c{-)/{p  +  fi),  while  the  discount  rate  is  f2/(/9  +  Jl). 

Optimal  stationary  policy  for  this  discrete-time  dynamic  programming  problem  is  gov¬ 
erned  by  the  Bellman  equations  [17], 


/?  + 


^  E  +  (n  -  v„  a  ,  j  =  0, 1, . . . ,  L  - 1 

-t- V'ijj,  i  =  l,2,...,  j  =  0,l,...,i-l 


(33) 

(34) 


where  Vjj  is  the  total  accumulated  cost  if  the  system  commences  with  i  customers  and 
with  the  interarrival  scale  being  j  [i.e.,  in  state  (t,  j)].  The  optimal  stationary  service  rate 
Hlj  for  state  (z,j)  is  then  the  minimizing  rate  for  the  corresponding  equation  of  (34).  If 
the  minimum  is  achieved  at  multiple  values  of  p,  the  lowest  rate  is  picked.  This  system 
of  equations  can  be  solved  numerically  via  the  value  iteration  method  [17].  Specifically, 
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beginning  with  Vij(O)  =  0  for  all  (i,  j),  we  iterate  the  system 


VoAk  + 1)  =  ^  E 

Vij{k  +  1)  =  ^min^  -^^1/1(0  +  c{n)  +  ^  ^  +  /^^t-ij(^) 

j'=i 


(35) 


(36) 


fib’) 

until  convergence.  The  optimizing  rate  at  each  iteration  k  then  converges  to  fi\y 

Based  on  this  approach,  optimal  server  control  is  designed  for  several  systems  with 
fractal  renewal  process  input,  and  the  resulting  policy  is  shown  in  Fig.  14.  To  allow  closed- 
form  minimization  of  (36)  in  each  step,  a  quadratic  service  rate  cost  c(/x)  =  is  adopted 
throughout.  Moreover,  a  linear  holding  cost  h(i)  =  hoi  is  used,  which  is  a  measmre  of 
expected  delay  in  a  first-in-first-out  queue.  In  the  cases  depicted,  cq  =  1  and  ho  =  0.01.  To 
reflect  equal  importance  of  the  present  and  the  future  in  om:  decision  making,  a  discount 
rate  0  =  0  is  used. 

Each  curve  in  Fig.  14  represents  the  optimal  service  rate  as  a  function  of  the  inter¬ 
arrival  scale  j,  with  the  number  of  customers  i  held  fixed.  For  a  given  queue  length,  the 
optimal  service  rate  is  seen  to  decrease  monotonically  for  coarser  (i.e.,  larger)  scales.  Thus, 
if  the  next  arrival  is  expected  to  be  distant,  the  service  can  be  relaxed  without  burdening 
the  future.  The  optimal  service  rate  is  also  monotonic  in  the  number  of  customers  when  the 
scale  is  fixed:  the  more  crowded  the  system,  the  busier  the  server.  In  Appendix  C,  we  show 
that  this  actually  holds  more  broadly  for  any  monotonic  convex  holding  cost  h{i).  Finally, 
we  note  that  small  7  cases  are  less  demanding  on  the  server.  In  effect,  heavier  tails  of  the 
corresponding  interarrival  densities  imply  that  prospective  customers  are  in  general  more 
distant,  and  hence,  have  less  impact  on  the  present  decision  making. 


5.1.2  State  Estimation  and  Realizable  Server  Control 


In  this  section,  we  develop  scale  estimates  that  can  be  used  together  with  queue  length 
information  in  the  state-based  controller  structure  of  the  previous  section  to  obtain  practical 
control  policies.  In  particular,  we  consider  estimation  of  the  interarrival  scale  from  the  time 
t  that  has  elapsed  since  the  last  arrival.  With  a  multiscale  representation  of  the  interarrival 
density  as  in  (2),  the  minimum  probability-of-error — i.e.,  maximum-arposteriori  (MAP) — 
scale  estimate  is  the  scale  j  that  maximizes 

^(j)  =  Pj  >t}  =Pj  exp  {-Xjt) . 


Using  (4),  we  get  that 

jC(j  -f-1)  //\  \\\ 

=  9  exp  (-(Aj+i  -  \j)t) , 

which  is  monotonically  increasing  in  t.  Thus,  C{j  -f  1)  >  C{j)  if  and  only  if 

where  the  equality  follows  from  (3).  Thus,  our  scale  estimate  j{t)  is  the  unique  integer  j 
such  that 


V 

(tj-  1)A’ 


i.e.. 


Kt)  = 


■  f  t{r)  -  l)Ay 


(37) 


Fig.  15  depicts  the  optimal  server  policy  of  Fig.  14(a),  modified  to  operate  without 
knowledge  of  the  interarrival  scale.  These  plots  are  obtained  by  warping  the  horizontal  axis 
to  yield  as  a  function  of  number  of  customers  present  i,  and  the  time  elapsed  since 
the  last  customer  arrival  t.  To  apply  the  prescription  of  Fig.  15,  we  follow  a  particular 
solid  curve  as  time  progresses,  until  a  customer  enters  or  leaves  the  system.  For  a  customer 
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arrival,  we  move  to  the  beginning  of  the  next  higher  curve.  Upon  a  customer  departxire,  we 
drop  vertically  to  the  next  lower  curve.  Thus,  past  history  of  the  system,  specifically  the 
arrival  epoch  of  the  last  customer,  is  crucial  in  our  decision  making. 

5.1.3  Simulations 

Through  simulations,  we  examine  the  benefits  of  our  multiscale  controller  over  more  conven¬ 
tional  designs  which  ignore  past  history,  such  as  the  M/M/1  controller  and  its  variations. 
In  contrast  to  our  systematic  multiscale  controller  design,  we  will  see  that  delicate  hand¬ 
tuning  is  generally  required  for  these  memoryless  controllers  to  achieve  even  reasonable 
performance  when  servicing  fractal  input. 

Our  simulations  are  based  on  the  next-event  time  advance  method  of  [12].  For  fair 
comparison,  each  controller  is  given  the  same  customer  input,  which  is  a  sample  function  of 
a  fractal  renewal  process  with  shape  parameter  7  =  1.8,  consisting  of  iV  =  500000  arrivals. 
The  cost  parameters  are  set  at  cq  =  10,  ho  =  0.01  and  =  0  throughout.  To  derive  the 
M/M/1  controller  of  [17],  the  input  is  treated  as  if  it  were  a  Poisson  stream,  with  arrival 
rate  estimated  according  to  N/Sx[N],  where  5a- [iV]  denotes  the  arrival  epoch  of  the  last 
customer.  The  resulting  service  policies  are  shown  as  the  dashed  plots  in  Fig.  16,  together 
with  the  multiscale  policies  which  are  represented  by  the  solid  curves. 

Based  on  the  simulated  systems,  we  have  obtained  the  total  cost  J  (defined  in  (30))  and 
the  average  delay  under  a  first-in-first-out  discipline,  and  have  tabulated  these  quantities  in 
Table  1.  From  the  table,  it  is  apparent  that  the  M/M/1  controller  is  substantially  inferior  in 
terms  of  either  measure.  As  we  can  see  in  Fig.  16,  this  controller  attempts  to  lower  overall 
cost  by  reducing  the  service  cost,  at  the  expense  of  a  higher  holding  cost.  Consequently, 
it  is  more  susceptible  to  congestion.  The  estimated  steady-state  customer  distributions  in 
Fig.  17  provide  further  evidence  for  a  longer  queue  expected  for  the  M/M/1  controller. 

One  way  to  lower  the  holding  cost  is  by  raising  service  rates.  This  can  be  achieved 
simply  by  scaling  up  the  service  rates  of  the  optimal  M/M/1  controller.  By  trial  and  error, 
we  have  tuned  the  M/M/1  controller  to  yield  the  least  overall  cost  for  the  given  customer 
input.  The  resulting  controller  is  depicted  in  Fig.  18,  alongside  the  optimal  multiscale 
controller,  and  its  performance  is  summarized  in  the  third  column  of  Table  1.  Although  it 


is  the  best  controller  obtainable  by  scaling  the  M/M/1  controller,  this  design  is  nevertheless 
still  suboptimal,  as  can  be  seen  from  Table  1.  Also,  as  is  apparent  from  Fig.  18,  the  modified 
M/M/1  controller  is  more  demanding  on  resources,  requiring  a  substantially  higher  peaJs 
service  rate  than  the  multiscale  controller. 

That  a  large  proportion  of  interarrivals  in  a  fractal  renewal  process  occur  on  short  time 
scales  suggests  that  an  alternative  controller  design  could,  in  principle,  be  based  on  mim¬ 
icking  fine-scale  behavior  of  the  optimal  multiscale  controller.  We  next  design  a  memoryless 
controller  of  this  type  based  on  the  fine-scale  service  rates  in  the  optimal  multiscale  policy. 
As  before,  we  hand-tune  these  service  rates  to  obtain  the  lowest  total  cost  for  the  given 
customer  input,  using  trial  and  error.  The  resulting  controller  is  plotted  in  Fig.  19,  while 
its  performance  is  summarized  in  the  last  column  of  Table  1.  As  is  apparent  firom  the  table, 
this  controller  is  still  inferior  to  the  multiscale  scheme;  note,  for  example,  the  considerably 
higher  average  delay. 

The  preceding  results  collectively  reflect  that  while  memoryless  service  policies  are  sim¬ 
ple  to  implement,  in  general  they  yield  rather  poor  results  for  fractal  customer  input.  Even 
when  caxefully  hand-tuned  according  to  the  actual  customer  arrivals,  these  policies  are  in¬ 
herently  inferior  to  the  multiscale  controller.  Because  the  fi-actal  model  is  a  process  with 
strong  memory,  intelligent  allocation  of  service  based  on  its  history  is  ultimately  key  to 
successful  control  of  this  class  of  queues. 

5.2  Flow  Control  for  Power-Law  Services 

In  this  section,  we  study  optimal  control  of  the  queueing  system  of  Section  4.4,  where 
Poisson  customers  are  serviced  with  power-law  holding  times.  The  basis  of  our  design  will 
be  the  transposed  multiscale  birth-and-death  process  of  Fig.  11,  which  is  used  for  queueing 
analysis  in  Section  4.4.  With  this  model,  it  is  the  service  that  is  described  in  the  form  of  a 
multiscale  representation. 

5.2.1  Multiscale  Flow  Control 

The  heavy-tailed  customer  distribution  inherent  in  such  queueing  systems  means  that  tra¬ 
ditional  flow  control  is  problematic.  In  terms  of  system  management,  such  controllers  result 


in  high  likelihood  of  buffer  overflow,  while  from  the  customers’  perspective,  long  delay  and 
poor  quality  of  service  are  experienced.  In  this  section,  we  describe  improved  flow  control 
policies  that  mitigate  these  problems.  To  achieve  this,  we  allow  the  controller  to  vary  the 
cidmission  rate  A  between  0  and  the  actual  arrival  rate  A.  As  before,  we  first  solve  this 
problem  assuming  complete  knowledge  of  the  state  of  the  system.  Thus,  we  denote  the 
input  rate  at  state  (i,  j)  by  Xij. 

The  objective  function  is  again  made  up  of  two  components,  with  a  holding  cost  h{i) 
which  reflects  buffer  occupancy,  and  a  throttling  cost  c(A)  which  penalizes  loss  of  customers. 
We  assume  that  the  holding  cost  is  monotonically  nondecreasing  in  the  number  of  customers, 
while  the  throttling  cost  is  monotonically  nonincreasing  in  the  admission  rate.  In  addition, 
we  assume  /i(0)  =  c(A)  =  0.  The  total  cost  is  accumulated  over  time,  with  a  discount  rate 
/3  included: 


(38) 

where  i{t)  is  the  number  of  customers  in  the  queueing  system  at  time  t,  and  X{t)  is  the 
admission  rate  at  time  t. 

Our  approach  to  this  problem  will  be  very  similar  to  the  server  design  of  Section  5.1.  Ap¬ 
plying  rate  uniformization  to  the  transposed  multiscale  birth-and-death  process  of  Fig.  11, 
we  obtain  the  equivalent  process  of  Fig.  20.  Note  that  here  we  have  also  lumped  the  zeroth 
superstate  into  a  single  state  to  obtain  a  more  realistic  model;  the  scale  of  the  next  service 
in  general  cannot  be  deduced  from  an  empty  queue.  We  next  recast  this  problem  into  its 
discrete- time  equivalent,  with  the  corresponding  Bellman  equations 


J  =  E 


’  poo 

/  e-^‘[h(z(t))-hc(A(i))]dt 


Vo  =  min 

A€[0,A]  P  + 


Vi  i  =  min  ^ 

A€[0,A]  + 


^  |c(A)  -h  A  +  (fl  -  A)  Foj  (39) 

—  {/i(i)  -I-  c(A)  -I-  ^  +  AVi+ij 

^  (40) 

+  z  =  l,2,...,  i  =  0,l,...,£-l. 


where  Vij  is  the  total  eiccumulated  cost  if  the  system  commences  in  state  (z,  j).  The  optimal 
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stationary  admission  rate  X*j  when  the  system  is  in  state  {i,j)  is  the  minimizing  rate  for 
the  corresponding  equation  of  (39)  and  (40).  If  the  minimum  is  achieved  at  multiple  values, 
we  pick  the  highest  rate.  In  contrast  to  the  optimal  server  control  problem  of  Section  5.1, 
the  zeroth-order  equation  (39)  is  no  longer  trivial,  since  the  throttling  decision  in  the  zeroth 
superstate  does  have  an  influence  on  the  future. 

Figs.  21  shows  the  optimal  stationary  flow  control  policy  obtained  by  solving  (39)  and 
(40)  via  the  value  iteration  method  [17].  Throughout,  a  quadratic  throttling  cost  and  linear 
holding  cost  are  used,  i.e.,  c(A)  =  co(A  —  A)^  and  h{i)  =  hoi,  where  cq  =  1  and  ho  =  0.01. 
Also,  the  discoimt  rate  is  set  at  (5  =  0. 

For  a  fixed  number  of  customers  i,  the  optimal  admission  rate  X*  j  is  seen  to  decrease 
for  coarser  scale  service  duration,  i.e.,  larger  j.  Thus,  if  the  current  job  is  expected  to 
require  long  service,  fewer  customers  should  be  admitted  to  prevent  high  costs  for  the 
futmre.  On  the  other  hand,  for  a  fixed  scale  j,  the  optimal  admission  rate  is  monotonic 
in  the  number  of  customers.  Hence,  higher  degree  of  throttling  is  required  for  a  busier 
system.  In  Appendix  D,  we  prove  that  this  monotonicity  actually  holds  more  generally  for 
any  convex  holding  cost  h{i).  Comparing  Figs.  21(a)  and  21(b),  we  also  see  that  lower  7 
cases  require  higher  throttling.  Due  to  the  longer  expected  service  time,  fewer  customers 
can  be  accommodated  in  these  cases.  Finally,  as  in  the  server  control  problem,  practicable 
flow  control  can  be  efiiciently  realized  from  the  idealized  policy  by  integrating  minimum 
probability-of-error  (MAP)  state  estimates: 


Kt)  = 


( <(^-  i)m\ 

\r}ln{l/q)J 


(41) 


where  t  now  represents  the  time  elapsed  since  the  last  service  completion. 


6  Conclusion 

In  this  paper,  we  have  applied  multiscale  techniques  in  the  study  of  fractal  point  processes  in 
various  realistic  networking  scenarios.  Our  analysis  of  these  processes  under  random  erasiue 
and  superposition  is  consistent  with  the  broadly-observed  self-similarity  in  aggregate  traffic; 
fluther  analysis  with  this  point  process  model  may  lead  to  further  insights  into  mechansims 
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by  which  self-similarity  arises  in  networks.  Our  queueing  analysis  resulted  in  methods  for 
computing  various  performance  measures  such  as  quality  of  service  and  resource  consump¬ 
tion.  This  analysis  identified  several  problems  with  the  use  of  traditional  queueing  design 
techniques  for  these  scenarios — in  particular,  the  substantial  underutilization  of  resources 
with  memoryless  service  of  fractal  point  process,  and  the  inherently  poor  quality  of  ser¬ 
vice  for  Poisson  process  serviced  with  power-law  holding  times.  To  mitigate  these  effects, 
we  have  applied  multiscale  paradigms  to  develop  systematic  design  methodologies  for  the 
operation  of  these  queueing  systems.  Several  practical  optimal  controllers  were  obtained 
from  this  approach.  Through  simulations,  we  have  shown  that  by  exploiting  system  history, 
a  simple  multiscale  controller  significantly  out-performs  simplistic  controllers  which  ignore 
such  past  information.  While  we  have  focused  on  several  specific  problems  mainly  for  illus¬ 
trative  purposes,  om  design  algorithms  are  readily  generalized  to  other  queueing  scenarios, 
such  as  those  involving  different  cost  structures. 


A  Derivation  of  the  Counting  Process  Distribution  Coeffi¬ 
cients:  Arrival-Observed  Case 

In  this  section,  we  derive  (t);  i  =  0, 1, . . .  | ,  the  arrival-observed  counting  process  distri¬ 
bution.  First,  as  given  in  (9),  a  closed-form  expression  exists  for  the  Oth-order  term 
To  obtain  Taylor’s  series  coefficients  for  higher-order  terms,  we  expand  the  ^-transform  of 
the  probability  distribution  (8)  as  follows 


p(2;  <)l’^  =  q  [  I  -I-  (-B  -I-  (At)  + 


(-B -h  2bTq)^  (At)2 
2! 


Extracting  the  coefiicient  of  z^,  we  get  that 

+  (42) 

where  for  /;  >  1, 

=  qB^-^b^ql"^  +  qB^'-^b'^qB!'^  4-  •  •  •  +  +  qb'^qB^'-^'^.  (43) 
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But  b  =  IB.  Thus, 

4^^  =  +  qB^^-^l'TqBlT  +  •  •  ■  +  qB^l^qB^-^l^  +  qBl^qB^'^^.  (44) 

This  is  simply  the  sum  of  all  terms  of  the  form 

qB™l'^qB”l'^ 


such  that  m  >  1,  n  >  0,  and  m  +  n  =  A:.  Moreover,  we  recognize  qB^l"*"  as  the  Zth  moment 
of  the  random  variable  R,  with  distribution  given  in  (12).  Thus,  we  have 

k 

=J2MiMk-i,  (45) 

J=i 


with  Ml  denoting  the  ith  moment  of  R. 

Next  we  turn  to  the  expansion  of  for  arbitrary  i.  First,  it  is  clear  that  the 

coefl&cients  of  are  zero  for  k  <i.  Thus,  we  have 


TT, 


(^)(t)  =  4)M 


i!  ‘^*+^(z  +  l)! 


(i) 


(At)*+i 


+ 


(46) 


Using  the  same  argument  as  before,  we  have  that  each  coefficient  is  the  sum  of  all  terms 
of  the  form 


qB"*i  iT  . . .  qB^^i  lTqB""‘+i  1^,  (47) 

with  mi,  m2, . . . ,  m^  >  1,  mj+i  >  0,  and  mi  +m2  -I - l-mj+i  =  k.  Such  terms  with  mi  =  i, 

where  1  <  Z  <  fc— Z  +  l,  are  those  with  m2, m3, . . .  ,mi  >  1,  mj+i  >  0,  m2H-m3H - l-mj+i  = 

k  — I .  But  these  are  exactly  the  terms  making  up  the  sum  ■  Thus,  we  have 

fc— i+l 

=  Y.  (48) 

i=\ 

B  Derivation  of  the  Counting  Process  Distribution  Coeffi¬ 
cients:  Random  Incidence  Case 

We  first  argue  that  the  left  null  spane  of  the  matrix  (— B  +  b^q)  is  spanned  by  the  vector 
qB“^.  Suppose  the  vector  v  is  in  the  left  null  space,  or 

V  (— B  +  b^q)  =  0 


where  0  is  a  zero  row  vector.  So, 


vB  =  vb^q  =  Acq 
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where  k  is  the  inner  product  vb^.  Right-multiplying  both  sides  by  B  \  then,  we  have 


V  =  K  (qB  . 


We  proceed  to  determine  the  coefficients  in  the  counting  process  probability  distribution 
i  =  0, 1, . . .  I .  Again,  as  given  in  (13),  we  have  a  closed-form  expression  for  the  0th- 

order  term  TrQ^(t).  For  higher-order  terms,  we  expand  the  z-transform  of  (8),  with  the  initial 
condition  assuming  the  value  p(z;  0)  =  a^qB"^.  Hence, 


p(z;  t)l"  =  d^qB  M  I  +  (-B  -f  zb™q)  (At)  + 


(-B -t- zbTq)2(At)2 
2! 


-1- 


1^. 


(49) 


As  before,  it  is  clear  that  the  coefficients  for  will  be  zero  for  k  <  i.  Thus,  we  have 


_  ji)  _  (i)  J _ 


(50) 


Focusing  on  the  A:th  power  of  t  in  the  coefficient  of  z*,  we  see  that  this  is  the  sum  of  all 
terms  of  the  form 


d^qB-^B’"!  lTqB""2  iT . . .  lTqB”‘-+i  1^,  (51) 

with  mi,  m2, . . . ,  mj  >  1,  mi^-i  >  0,  and  mi  +  m2  -!-•••-)-  mi+i  =  k.  For  those  terms  with 
mi  =  1,  the  product  reduces  to 

d2qBm2  iT  , . .  qgmi  iTqjgmi+i  ^52) 

since  qlT  =  1.  But  the  sum  of  all  terms  of  this  form  is  just  d^o^*  On  the  other  hand, 
terms  with  mi  >  1  can  be  written  as 


dSqB^i  lTqB”»2  iT  . . .  qB"**  lTqB"*i+i  1^,  (53) 


with  mi,  m2, . ..  ,mi  >  1,  mi+i  >  0,  and  tui  -|-  m2  -I - h  mj+i  =  A;  —  1.  But  this  is  precisely 

d^a^^li-  Thus,  we  have  shown  the  relation 

’’fc  ^  =  (4-1  +  (54) 

C  Optimal  Fractal  Queue  Server  Monotonicity 

In  this  appendix,  we  show  that  the  optimal  service  rate  fx*j  for  a  fractal  renewal  process  input 
satisfies 


for  i  =  1,2, ... , 


(55) 


if  the  holding  cost  function  h{i)  is  convex,  i.e., 

h{i)  —  h{i  —  1)  <  h{i  4- 1)  —  h{i),  for  i  =  1, 2, - 

We  first  show  that  (55)  holds  if  the  first  difference  of  Vij,  defined  as 

^i,j  —  ^i,j  l,j) 

is  nondecreasing  in  i.  Rewriting  the  Bellman  equations  (34)  for  i  >  0,  we  get  that 


(56) 


+  mm 
#'€[0,71] 


cifi)  -  /iAi 


1. 


Now,  suppose  that  Ajj  is  indeed  nondecreasing  in  i,  and  let  fi,i2  be  some  nonnegative  integers, 
with  i2>  h-  Then,  for  every  ju  < 


<  c{n)  —  /lAij  j  j(At2,j  ~  Aij_j) 

<  c{fi)  —  —  Aij  j)  =  c{fi)  — 

where  the  first  inequality  follows  from  the  definition  of  //*j  j,  and  the  second  inequality  follows  from 
monotonicity  of  Ajj.  So,  j 

We  proceed  to  show  that  Aij  is  indeed  monotonic  in  i.  For  this,  it  suffices  to  show  that  at  each 
stage  k  of  the  value  iteration  method,  the  first  difference 

Aij(fc)  =  yij(fc)-Vi-i,i(fe)  (57) 


has  this  property.  We  set  the  boundary  value  AQj{k)  to  0,  for  all  j,k. 

Now,  it  is  trivial  that  Aij(O)  is  nondecreasing  in  i.  Assuming  this  holds  for  Aij(k),  we  consider 
the  next  iteration.  Specifically,  we  rewrite  the  value  iteration  equation  (36)  as 


m 


(58) 


+  min 


c(fi)  -  fiAi,j(k) 


,  i  —  1, 2, ... . 


Defining  the  optimal  service  rate  at  each  stage  k  of  the  value  iteration  method  as  fijj,  i.e.. 


(k)  A 

fz}  /  =  arg  min 


c(fx)  -  nAi,j{k) 


(59) 
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for  every  i,j,  k,  we  have  that 


Ai+ijik  +  1)  =  Vi+ij{k  +  1)  -  Vij{k  +  1) 

-  +  +  :;^J2^^V-^Vi+2,f{k) 

-  Hi)  -  cr^q^''^Vi+ij>{k) 

^  i'=i 

-  («  -  ^)  y^Ak)  - 

=  ^  E  AA-^Ai+,,j.ik) 

+  (n-  {Ai+i,j{k)  -  Ai,,(fc))|. 

Similarly,  we  get  that 

Ai,jik  +  1)  <  -  Ki  -  1)  +  S  cr^qi'-^Ai+ij'ik) 

+  (n-  Aijik)  -  (Aijik)  -  A,_x.^(A:))|. 

Subtracting  these  two  inequalities,  we  get  that 

{/3  +  Q){Ai+ij{k  +  1)  -  Aij{k  +  1))  >  ^{h{i  +  1)  —  h{i))  -  {h{i)  -  h{i  —  1)) 

A  ^ 

■'■  'iiFi  X)  Aqi'-'^\Ai+2,j-  (A:)  -  Ai+ij-(A:)] 

^  i'=i 

~  7ji-l  ~  ['^i+l,j(^)  ~  '^i,j(A^)] 

+  (*)]^  • 


(60) 


(61) 


(62) 


By  the  induction  hypothesis,  convexity  of  h(i),  and  the  fact  that  fl  >  X/if~^  -  Ai+i  j,  we  get  that 
Ai+ij{k  +  1)  >  Aij{k  +  1).  Our  assertion  therefore  follows  by  mathematical  induction. 


D  Optimal  Flow  Control  Policy  Monotonicity 

In  this  appendix,  we  show  that  the  optimal  admission  rate  A*j  for  a  power-law  queueing  system 
servicing  a  Poisson  input  satisfies 

fori  =  1,2,...,  (63) 
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given  that  the  holding  cost  function  h{i)  is  convex,  i.e., 

h{i)  —  h{i  —  1)  <  h{i  +  1)  —  h{i),  for  i  =  1, 2, - 

Using  exactly  the  same  argument  as  in  Appendix  C,  we  can  show  that  A*  j  is  nonincreasing  in 
i,  if  the  first  difference  of  Vij  is  nondecreasing  in  i.  We  proceed  to  show  that  the  first  difference 
Aij,  defined  as  in  Appendix  C,  is  monotonically  nondecreasing.  Our  approach  will  exploit  the  value 
iteration  method,  and  prove  by  mathematical  induction  that  the  first  difference  of  Vij  (k)  at  each 
stage  of  the  iteration,  is  monotonic  in  i.  Again,  Ao{k)  is  set  to  be  zero  for  all  k.  Also,  it  is  dear  that 
Ai,j(0)  is  nondecreasing  in  i,  since  V'ij(O)  is  identically  zero.  We  assume  the  assertion  holds  for  k, 
and  analyze 


+  1)  —  Ai_ij(fc  +  1). 

=  Voik),  and  Aoj(fc)  =  . 

To  this  end,  we  write  the  value  iteration  equations  as 


For  each  k,  we  also  define  =  Vo{k),  and  Ao,j{k)  =  Ao(fc). 


Vb(A:  +  l)  = 


^  ^^oik)  +  mill, 
P  +  “  I  A€[0,A] 


c(A)  +  A  aV'~'  (ViAk)  -  Voik)) 


j'=i 


,(fc) 


(64) 


+  min_  [c(A)  +  AAi+i,j(fc)]l,  i  =  l,2,...,  j  =  0, 1, . . . ,  L  -  1. 

A6[0,A]  J 

We  define  the  optimal  service  rate  at  each  stage  k  of  the  value  iteration  method  as  i.e.. 


(65) 


.(Ic)  A 

A>  :  =  arg  mm 

A€[0,A]  L 


c(A)  +  AAi+i,j(A:) 


(66) 


Now,  for  i  >  1,  we  get  that 


and 


V5«.,-(l=  + 1)  =  + 1)  +  +  (tl  -  ^) 


'•  '  j'=l 

+  c(A«)  +  Al5A,+i,^(fc)|. 


»+i>j 


(k) 
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Subtracting,  we  get  that 


Aj+i j(A;  +  1)  =  Fi+i,_,(/:  +  1)  -  Vij{k  +  1) 

+  (fi  -  Ai+iAk)  +  a1+\.,-  (Ai+2j(fe)  -  Ai+i,,(/:))|. 

Similarly, 

Aijik  +  1)  =  ViAk  + 1)  -  Vi.^Ak  + 1) 

-  ;9Tn{'‘<')-'‘<‘'-i)  +  sfe'E‘’V'-‘A,-.j.(4)+  (n-^)  AyW 

+  Al?,j(AM.u(i=)- A, „■(*))}■ 

Thus, 

+  n){Ai+iAk  + 1)  -  AiAk  + 1))  >  (^{hii  + 1)  -  hit))  -  {h{i)  -  Hi  - 1)) 

J-i  (67) 

+  ^  [Ai+i.,(fe)  -  A,,,(A;)] 

+  AW,,.[A,+2.,(fe)-Ai+i,,(A:)]). 

Thus,  it  follows  from  the  induction  hypothesis,  convexity  of  h(i),  and  the  fact  that  Q  >  + a|*\  ^  ^ 

that  Aij{k  +  1)  is  nondecreasing  for  z  >  2.  To  complete  the  proof,  we  now  show  that  A2,j{k  +  1)  > 
Aij{k  +  1).  We  note  that 

A2Ak  +  1)  =  V2Ak  +  1)  -  ViAk  +  1) 

^  +  5^  E'V'-a.,,,w  +  (fi-  A,,,(t) 

+  4'J(A,,,(fc)-A,^W)|, 

and 

AiAk  +  1)  =  ViAk  +  1)  -  Vo,j{k  +  1) 

+  Ag(A2,^(fc)-Ai,,(A:))|. 
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Thus,  again,  our  assertion  follows  for  the  case  f  =  1,  for  the  same  reasons. 
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Figure  1:  Successive  magnification  of  the  counting  process  associated  with  a  fractal  renewal 
process;  (a)  the  original  process;  (b)  zoomed  version  of  (a);  (c)  zoomed  version  of  (b). 
Interarrivals  were  synthesized  according  to  the  power-law  density  function  (1),  with  shape 
parameter  7  =  1.5 
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0  arrivals  1  arrival  2  arrivals 

Figure  2:  The  multiscale  pure-birtb  process  based  on  a  hnite-scale  representation;  dashed 
boxes  denote  conceptual  grouping  into  superstates.  The  rate  of  departure  from  each  state  is 
a  function  the  scale,  Xj  =  X/rf;  the  probability  of  entering  scale  f  of  succeeding  superstate 

IS  pji  =  (7  . 


Table  1:  Performance  of  various  queueing  server  controllers  servicing  fractal  renewal 

process  input 


Multiscale 

M/M/1 

Modified  M/M/1 

Modified  Fine-Scale 

Controller 

Controller 

Controller 

Controller 

Overall  Cost 

0.9963e6 

1.161e6 

1.009e6 

1.017e6 

Average  Delay 

44.4 

115.0 

46.8 

49.0 
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Figure  3:  First  4  terms  in  the  arrival-observed  counting  process  probability  distribution 
for  a  fractal  renewal  process,  computed  with  a  20-scale  dyadic  representation;  the  shape 
parameter  of  the  interarrival  distribution  is  7  =  1.8.  The  time  axis  is  normalized  with 
respect  to  the  hnest-scale  arrival  rate  X. 


Figure  4:  First  4  terras  in  the  random  incidence  counting  process  probability  distribution 
for  a  fractal  renewal  process,  computed  with  a  20-scale  dyadic  representation;  the  shape 
parameter  of  the  interarrival  distribution  is  j  =  1.8.  The  time  axis  is  normalized  with 
respect  to  the  finest-scale  arrival  rate  X. 


Figure  5:  Interarrival  density  function  of  a  fractal  renewal  process  under  Bernoulli  era¬ 
sure,  with  erasure  probability  p.  The  shape  parameter  of  the  original  process  is  7  =  1.8. 
These  computations  were  performed  with  a  20-scale  dyadic  representation.  The  plots  are 
normalized  with  respect  to  the  hnest-scale  arrival  rate  A. 
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Figure  6:  Arrival-observed  counting  process  distribution  of  the  superposition  of  two  in¬ 
dependent  fractal  renewal  processes,  generated  with  the  counting  process  results  of  Figs.  3 
and  4.  The  shape  parameter  of  both  constituents  is  7  =  1.8.  The  time  axis  is  normalized 
with  respect  to  the  hnest-scale  arrival  rate  A. 


Figure  7:  Random  incidence  counting  process  distribution  of  the  superposition  of  two 
independent  fractal  renewal  processes,  generated  with  the  counting  process  results  of  Fig.  4. 
The  shape  parameter  of  both  constituents  is  7  =  1.8.  The  time  axis  is  normalized  with 
respect  to  the  Bnest-scale  arrival  rate  X. 


0  in  system  1  in  system  2  in  system 

Figure  8;  The  multiscale  birtb-and-deatb  process  based  on  a  bnite-scale  representation. 
This  is  obtained  by  adding  death  transitions  to  the  multiscale  pure-birtb  process  of  Fig.  2. 
To  model  a  single-server  memoryless  queueing  system,  all  death  transitions  occur  at  equal 
rate  p. 
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Steady-state  probabilities 


Figure  9;  First  21  terms  in  the  steady-state  customer  distribution  corresponding  to  a 
memoryless  single-server  queueing  system  with  fractal  customer  arrivals.  The  input  process 
is  modeled  with  a  20-scale  dyadic  representation,  with  7  =  1.8,  and  Ijp  =  p/\  =  0.8. 
The  solid  curves  represent  simulation  results  obtained  with  the  next-event  time  advance 
simulation  [12],  while  the  dashed  curves  represent  computed  results. 


43 


Decay  rate, 


Figure  10:  Decay  rate  of  the  steady-state  customer  distribution  in  a  single-server  queueing 
system  with  fractal  renewal  process  input.  The  solid  curves  are  obtained  via  bisection- 
search  solution  to  (26).  The  dashed  curves  are  the  closed- form  upper  bound  prescribed  by 
(28).  A  20-scale  dyadic  representation  is  used  in  this  computation. 
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0  in  system  1  in  system  2  in  system 

Figure  11:  The  transposed  multiscale  birtb-and-deatb  process  for  memoryless  input  ser¬ 
viced  by  power-law  server.  This  process  is  obtained  by  reversal  of  tbe  roles  of  birtbs  and 
deaths  in  tbe  process  of  Fig.  8. 
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Steady-state  Probabilities 


Figure  12:  Comparison  of  simulated  and  theoretical  steady-state  customer  distribution  for  a 
queueing  system  with  power-law  bolding  time  servicing  memoryless  input.  The  solid  curves 
represent  results  of  discrete-event  simulation  of  the  queueing  system.  The  dashed  curves 
represent  theoretically-predicted  distributions  for  the  corresponding  queueing  scenarios.  For 
the  case  7  =  1.5,  p  =  X/p,  =  2  x  10“^,  while  for  the  case  7  =  1.8,  p  =  X/p  =  1  x  10“®. 
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£2— Xfl  £2— 2.g  — |j,jQ  £2  2,g  — (12,0 


^  ^t-l  ^  ^L-l  ^  ^Z,-l  ^*'2,^.-l 

0  in  system  1  in  system  2  in  system 

Figure  13:  Tie  continuous-time  Markov  process  employed  in  server  controller  design.  This 
process  is  obtained  from  the  multiscale  birtb-and-death  process  of  Fig.  8  by  adding  self¬ 
transitions  such  that  the  total  rate  leaving  each  state  is  ft.  Note  that  system  dynamics  are 
preserved. 
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Figure  14:  State-based  multiscale  service  policy  for  a  queueing  system  with  fractal  customer 
arrivals.  The  shape  parameter  of  the  input  is  7  =  1.8  in  (a)  and  7  =  1.2  in  (b).  The  holding 
cost  h{i)  =  O.Oli,  service  rate  cost  c(p)  =  and  discount  rate  /3  =  0  are  used.  Note  that 
with  no  customer  present,  optimal  service  rate  is  identically  0. 


48 


0.5 


Figure  15:  Realizable  optimal  server  for  fractal  renewal  process  input,  based  on  the  number 
of  customers  in  the  system,  and  the  time  elapsed  since  the  last  arrival.  The  depicted  policy 
is  for  the  case  7  =  1.8,  c(ju)  =  p?,h{i)  =  O.Oli,  and  is  obtained  via  a  rewarping  of  the 
horizontal  axis  in  the  graph  in  Fig.  14(a). 
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Figure  16:  Stationary  server  control  policies  for  queueing  systems  servicing  fractal  renewal 
process  input,  for  the  case  7  =  1.8,  c(p)  =  lOp^,  h{i)  =  O.Oli,  /?  =  0.0.  The  solid  curves  rep¬ 
resent  the  optimal  multiscale  server.  The  dashed  curves  denote  an  M/M/1  service  controller 
designed  with  the  input  treated  as  a  Poisson  process. 
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steady-state 


r  o  mipiieing  system  arising 

c(^)  =  lOpl^  h{i)  =  P 


0.4 


Figure  18:  Stationary  server  control  policies  for  queueing  systems  servicing  fractal  renewal 
process  input,  for  the  case  7  =  1.8,  c{p)  =  lOp^,  h{i)  =  O.OH,  /3  =  0.0.  The  solid  curves  rep¬ 
resent  the  optimal  multiscale  server.  The  dotted  curves  denote  the  policy  for  the  optimized 
variation  of  the  M/M/1  queueing  controller  obtained  by  scaling  the  service  rates. 


-  Multiscale  controller 

-  Modified  fine-scale  controller 


Figure  19:  Stationary  server  control  policies  for  queueing  systems  servicing  fractal  renewal 
process  input,  for  the  case  7  =  1.8,  c(p)  =  lOp^,  h{i)  =  O.Oli,  S  =  0-0.  The  solid  curves 
represent  the  optimal  multiscale  server.  The  dot-dashed  curves  denote  the  optimized  policy 
obtained  by  scaling  the  hne-scale  service  rates  taken  from  the  multiscale  policy. 


0  in  system 


^  M-L-1  ^1,L-1  ^  ^  2,i-l 

1  in  system  2  in  system 


Figure  20:  The  continuous-time  Markov  process  used  in  how  control  policy  design.  This 
process  is  obtained  from  the  transposed  multiscale  birth-and-death  process  of  Fig.  11  by 
adding  self-transitions  such  that  the  total  rate  leaving  each  state  is  Q.  Also,  the  zeroth 
superstate  is  lumped  into  a  single  state. 
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Figure  21:  State-based  multiscale  Slow  control  for  a  queueing  system  with  power-law 
service  time  and  Poisson  customer  arrivals.  The  shape  parameter  of  the  service  duration 
distribu^on  isj  =  1.8  in  (a)  and  7  =  1.2  in  (b).  Holding  cost  h{i)  =  O.Oli,  throttling  cost 
c(A)  =  (A  —  A)^,  and  discount  rate  jd  =  0  are  used. 


